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Fast Molecular-Dynamics Simulation for Ferroelectric Thin-Film Capacitors Using a 
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A newly developed fast molecular-dynamics method is applied to BaTiOs ferroelectric thin-film 
capacitors with short-circuited electrodes or under applied voltage. The molecular-dynamics simu- 
lations based on a first-principles effective Hamiltonian clarify that dead layers (or passive layers) 
between ferroelectrics and electrodes markedly affect the properties of capacitors, and predict that 
the system is unable to hop between a uniformly polarized ferroelectric structure and a striped ferro- 
electric domain structure at low temperatures. Simulations of hysteresis loops of thin-film capacitors 
are also performed, and their dependence on film thickness, epitaxial constraints, and electrodes are 
discussed. 
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I. INTRODUCTION 

Ferroelectric thin films are beginning to see wide- 
ranging applications, for example in multilayer capaci- 
tors, nonvolatile FeRAMsji and nanoactuators. There is 
strong pressure to reduce the sizes of such thin-film struc- 
tures. In recent years, the preparation of oxide thin films 
by low-temperature non-equilibrium techniques such as 
molecular beam epitaxy and pulsed-laser deposition have 
attracted a great deal of attention, as they enable finely 
controlled growth of epitaxial thin films. '^^ 

It is well known that the properties of ferroelectric ca- 
pacitors are highly influenced by the properties of the 
interface between the ferroelectrics and the electrodes. 
For example, the fatigue of ferroelectric capacitors is as- 
sociated with the appearance of dead layers (or passive 
layers) near the electrodes^'i^i^i^ and imperfect electrodes 
cannot fully screen the polarization of ferroelectrics;^ 
leading to a finite depolarization field in the ferroelectric 
film. However, the nanosize effects and temperature de- 
pendences of ferroelectric capacitor hysteresis, polariza- 
tion switching, and dynamics of domain wall motion re- 
main poorly known. Experimentally, in situ observations 
are difficult. Theoretically, the long-range Coulomb in- 
teraction limits the size and time of molecular-dynamics 
(MD) simulations, and it has been unclear how to in- 
clude surface effects and depolarization fields caused by 
interface structures. 

In 1994, King-Smith and Vanderbilt studied the total- 
energy surface for zone-center distortions of perovskite- 
type ferroelectric oxides ABO^ {A is a monovalent or 
divalent cation and i? is a penta- or tetravalent metal) at 
zero temperature using first-principles calculations with 
ultrasoft-pseudopotentials and a plane- wave basis seti^ 
Starting from the full symmetric cubic perovskite struc- 
ture, they define the displacements w^ of atoms r {—A, 
B, Oi, On, Oiii) in the Cartesian directions a{= x,y,z) 
along the Fis soft-mode normalized direction vectors ^a 



\v2' 



^QSa ~~ ^a 



/ 



Sa 
;^Oii 



(1) 



with the scalar soft-mode amplitude Uq. Under the con- 
dition that the strain components rji (i — 1, . . . , 6; Voigt 
notation; rji ~ en, 774 — 623) minimize the total energy 
for each u — {u^, Uy,Uz), they expressed the total energy 
as 

^tot ^ £;0 ^ ^^2 ^ ^,^4 ^ y (y2 2 ^ 2^2 _^ ^2^2) ^ (3) 



where u^ 



\, E^ is the total energy for 



the cubic structure, k is half the eigenvalue of the soft 
mode, and a' and 7' are the constants determined from 
coupling constants between atomic displacements and 
strains. Their expression properly describes the coupling 
of polar atomic-displacement and strain degrees of free- 
dom. 

In 1994-1997, Zhong, Vanderbih, and Rabei^iii and 
Waghmare and Rabe^^ expanded Eq. ^ from a mean- 
field framework to a local-mode framework, replacing u 
by {u\, where the braces {} denote a set of it in a simu- 
lation supercell, as 

+ V°'^^(r?i,---,%)-Hr"*(M,ryi,...,,76) ■ (3) 

Here y^^'^, y^P', y'''^°'*, y°^^^ and y'"' are a local- 
mode self-energy, a long-range dipole-dipole interaction, 
a short-range interaction between soft-modes, an elas- 
tic energy, and an interaction between the local modes 
and local strain, respectively. They employed Eq. ([3]) 
as an effective Hamiltonian for {u\ in the supercell. 



performed Monte-Carlo simulations, and demonstrated 
the ability to describe the phase transitions of bulk 
ferroelectrics. The coarse- graining that reduces the 
15-dimensional atomic displacement vector v^^ to a 3- 
dimensional local soft-mode amplitude vector Ua in each 
unit cell was shown to be a good approximation. How- 
ever, the computation of V^^^ was still time-consuming, 
owing to the long-range Coulomb interaction, thus lim- 
iting system size and simulation time that could be han- 
dled in practical simulations. 



In 2003, Waghmare, Cockayne, and Burton introduced 
a technique to decrease the computational time for y^P' 
(or forces exerted on {u})J^ Direct calculation of the 
forces in real space requires a computational time pro- 
portional to N'^, i.e., 0(N'^), where TV is the supercell 
size (N = L.J. X Ly X L^). It decreases to O(iVlogiV) if 
one calculates the forces in reciprocal space using fast- 
Fourier transform (FFT) methods. This acceleration in 
computational speed enabled us to perform MD simula- 
tions on {u} in a large supercell, and was applied to bulk 
relaxor ferroelectricsj^^ii^ 



Here, we explain how the fast MD method for simulat- 
ing a first-principles effective Hamiltonian can be applied 
to study ferroelectric thin-film capacitor structures with 
short-circuited electrodes or external electric fields. This 
new MD method can simulate perovskite-type ferroelec- 
tric thin-film capacitors with dead layers and consequent 
depolarization fields. The high speed of this MD method 
enables us to simulate a ferroelectric material for a real- 
istic system size (up to 100 nm) and a realistic time span 
(> 1 ns). 

In the next section, we explain the formalism of the 
new MD-simulation method for thin-film capacitors. Re- 
sults of simulations of BaTiOa bulk and thin-film capaci- 
tors are shown in Sec. lIIII In subsection llll Al we confirm 
the reliability of our MD program by simulating thermal 
properties of bulk BaTiOa. The advantage of this MD 
method compared to the Monte-Carlo method is also dis- 
cussed. In subsection IIIIBi we perform heating-up and 
cooling-down simulations for thin-film BaTiOs capacitors 
with perfect and imperfect electrodes. Thickness depen- 
dence of simulated striped domain structures in thin-film 
capacitors with imperfect electrodes are analyzed in de- 
tail. We have already reported some simulated results 
of thin-film capacitors of this subsection and determined 
thermal properties in Ref. [l5| briefly. In subsection llll C[ 
newly obtained simulated results of hysteresis loops of 
thin-film capacitors are reported. In Sec. IIV[ we summa- 
rize the paper. 
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II. FORMALISM AND METHOD OF 
CALCULATION 

A. Effective Hamiltonian 

The effective Hamiltonian used in the present MD sim- 
ulations is basically the same as that in Ref. [i3|. Here, 
we present the Hamiltonian with a notation similar to 
that in Ref. [IM as 



H"" ^ 



M* 



2 ^-^ 

R,a 



uliR) 



M. 



acoustic \ ^ ^;,2 

2 ^ 

R.a 



^liR) 



+ ^^™"P' ''°'"°({m}, m,--;Ve) + ■^™"P' '"''°({m}, {w}) 

^z*Y.e-u{R), (4) 



R 



where u = u{R) and w ~ w{R) are, respectively, the 
local soft-mode amplitude vector and the local acoustic 
displacement vector of the unit cell at i?, the a compo- 
nent of R runs over 



Ra — 0, ao, 2ao, 



{La - l)ao , 



(5) 



M. 



■,riQ are the homogeneous strain components, and 
and M* „„„Hn are the effective masses for u and 



dipolc ^^^"-^ ^''-^acoustic 

w, respectively. Note that u can also be considered as 
the optical displacement, in contrast to the acoustic dis- 
placement ly, or the dipole moment Z*m, where Z* is the 
Born effective charge associated with the soft mode. In 
the effective Hamiltonian ((H), external electric field £ is 
taken into account through its vector product with each 
dipole moment Z*u. 

To determine the effective mass M^jp^j^,, let el^{k,i) be 
a mass-weighted i-th eigenvector of the phonon dynam- 
ical matrix at wavevector k. Its eigenvalue {a;(fe,z)}^ 
is the corresponding phonon frequency. Moreover, let 
(ijj,(fc, i) — e^(fc, i)/\/MT be an atomic displacement vec- 
tor, which is normalized as X^ari'^aC'''*)}^ ~ ^ ^y ^'^" 
justing the norm of e{k,i). Here, Af^ is the mass of 
atom r. Generally, the effective mass of a phonon is k- 
and mode-dependent: 



M*{k, 



\dl{k,^)}H^ 



(6) 



However, as an approximation, we have to employ a 
unique effective mass for dipoles in the MD simulation. 
Thus using the steepest descent Fis soft-mode normalized 
direction vectors 1^ = (0.20, 0.76,-0.21,-0.21,-0.53) 
and ^x = ^y = from Ref. [ll|, for BaTiOa, we set 
^^dipoic as 



M, 



ihttp://lQtQ . SQurcef orge.net /fercun/y 



dipole 



\CzVAr = 39.0 amii 



(7) 



It should be mentioned that ^J is not equal to the d^ 
of the Fis soft-mode of phonon, because M^, M^ , and 
M^ are not identical. 

The local-mode self-energy V^'^^^{{u}) is 



N 



7 [ul{R,)ul{Ri) + ul{Ri)ul{Ri) -t- ul{R,)ul{R,)\ } , 

(8) 

where u^{Ri) = ul{R,) + ul{Ri) + ul{Ri). 

The long-range dipole-dipole interaction V^'^p'({m}) is 

-. N N 

V^'P'(M) ^ ^Y.T.T.T.''o.{R^)<^o.p{R^,)u,{R,) , 



i=l a j=l p 



where 



(9) 



^MR.^) - — L iR-TTA^ ' 

(10) 
Eoo is the optical dielectric constant (or refractive index 
squared), dap is the Kronecker delta, a hat indicates a 
unit vector, n is the supercell lattice vector 



— 2LaaQ, —LattQ, 0, Latto, 2Laao, 



and flo is the equilibrium lattice constant. In Eq. (jlOp . 
y^ indicates that the summation does not include terms 
for which Rij — n — 0. 

We take account of short-range interactions between 
the optical displacements u(R) up to third nearest neigh- 
bor (3nn) as 

N 3nn 

y^'^-*(M) = -^^^^«„(ilO J,,,„/3^/3(il,) , 

i—1 a j (3 

(12) 
where Jij^ap is the short-range interaction matrix, which 
can be classified into 7 independent interaction parame- 
ters;^ Jtj^ap = ±jk {k = !,■■■ ,7). 

In practice, K2uf in Eq. ^, Eq. (O, and Eq. (fT2l) . 
in which Ua is quadratic, are gathered and calculated in 
reciprocal space as 

v^^^^iiu}) ^ i^^^;(fc)$r;'(fc)"/3(fc), (13) 

k a,p 

where Ua{k) is the Fourier transform 

Ua{k) = y^ Ua{R) exp(— ifc • R) , 

R 



(14) 



of Ua{R), ^'^p (fc) is similarly the Fourier transform of 
the quadratic interaction matrix (which is only calculated 



once at the beginning of simulation), ^'^ and A; is a recipro- 
cal vector in the first Brillouin zone of the unit cell such 



fha 



La~12TT 



1 27r 1 27r 1 2tt 
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Lot CiQ ^a ^0 ^ ^0 

(15) 
The homogeneous elastic energy v°^^^'^°'^°{rii, ■ ■ -,776) 



2La flo 



IS 



V 



clas, homo 



{Vir--,V6) = y^ii(?7?+'72 +^3) 



+ NBi2{rj2r]3 + 773771 + 771772) 
+ jB^iirjl+Tjl+Tjl) , (16) 

where -Bii, -B12, and B44 are the elastic constants ex- 
pressed in energy unit {Bn = QqCu, B12 — aQCi2, and 
-B44 — 09(744). 

The inhomogeneous elastic energy V°^^'^'^^°{{w}) is 
also calculated in reciprocal space as 

fc a,P 

(17) 
For the force constant matrix <I*° «^''" °{k), we employed 
the long-wavelength approximation. For instance, the 
diagonal part is 



$' 



;las, inlio 



1 



(fc) = M [^-^11 + ^a^44 + fc,'B44] 



N 



(18) 



(11) and the off-diagonal part is 



$ 



clas, inho 
xy 



(fe) = — [k^kyBi2 + k^kyBii] . (19) 



The coupling between {u} and homogeneous strain is 
the same as that given in Ref. Q, i.e., 

^coup,ho.no(|^|^ m, •••,%) = i ^ ^ ^ 77. a, y,{R) . 



R i = l j=l 



(20) 

Here, yi(i2) = ul{R), V2{R) = ul{R), y^iR) = ul{R), 
Va{R) = Uy{R)uz{R), ysiR) = Uz{R)ux{R), and yeiR) 

= Ux{R)Uy{R), 



C = 



(21) 

and Bixx, Biyy, and B^y^ are the coupling coefficients 
defined in Ref. [3]. 

The coupling between {u} and inhomogeneous strain 
is also calculated in reciprocal space as 

k a i—1 

(22) 
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where Wa{k) and yi{k) are the Fourier transforms of 
Wa{R) and yi{R), respectively. For the 3x6 couphng 

I 



matrix B(fe), we again employed the long-wavelength ap- 
proximation 



B(fe) ^ ^ 



^x^lxx '^X^lllll f^x^ 
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(23) 



In the present MD simulations of BaTiOa, the param- 
eters from Refs. [ifll and [il|, which are determined by 
first-principles calculations, are employed. As mentioned 
in Refs. [10| and [ll|, this parameter set leads to an un- 
derestimation of the Curie temperature Tq. To correct 
this underestimation, we follow these references in apply- 
ing a negative pressure of p = —5.0 GPa in all simula- 
tions. 



B. Molecular Dynamics 

MD simulations with the effective Hamiltonian of 
Eq. ([4]) are performed in the canonical ensemble using the 
Nose-Poincare thermostat. ^^ This simplectic thermostat 
is so efficient that we can set the time step to Af = 2 fs. 
In our present simulations, we thermalize the system for 
40,000 time steps, after which we average the properties 
for 10,000 time steps. 

In Fig. [1] we roughly illustrate how to calculate the 
forces exerted on Ua{R) with <&^T (fc) in Eq. (fT3|) and 
how the time evolution is simulated. First, Ua{R) is 
FFTed to u^ik), the force F„(fc) = - E /3 ^Tp'^{k)ufi{k) 
is calculated in reciprocal space, and then the force in 
real space is obtained by the inverse FFT of Fa{k). In 
practice, updates of Ua{R) and Ua{R) — ^Ua{R) are 
processed in the manner of the Nose-Poincare thermo- 
stat. 

The homogeneous strain components ?7i, • • •, jye are de- 
termined by solving 






clas. homo / 



+ ^^°"P''^°'"°(M,m, •••,%)] =0 (24) 



at each time step according to {u\ so that »7i, • • •, Jye niin- 

imize l/<=las,homo(^^^ . . .^^g)_,_^coup,homo(-|^| ,^^^ . . .^^^^ 

While the local acoustic displacement Wa{R) could be 
treated as dynamical variables using the effective mass 
-^acoustic: ^^ havc instead chosen to integrate out these 
variables in a manner similar to the treatment of the 
homogeneous strain. That is, Wa{R) is determined 
so that yeias,inho(|^|) _^ v"™"?^ '"''°({m}, {lo}) bccomcs 
minimum at each time step according to Ua{R). Techni- 
cally, the minimization is performed by solving the linear 
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FIG. 1: Simplified fiow chart for calculating forces on Ua{R). 
Fast Fourier transform (FFT) and inverse FFT (IFFT) enable 
rapid calculation of long-range dipole-dipole interactions. 



set of equations 

Jelas, inho(^^ ~(;^) ^ B(fc)y(fe) = (25) 

for each k in reciprocal space. 

C. Ferroelectric Thin Films 

If a ferroelectric thin film is placed in isolation in vac- 
uum without electrodes as depicted in Fig.[2l^a), its spon- 
taneous polarization P = {Px,Py,Pz), which is repre- 
sented by a thick arrow in the figure, induces charges 
icTjnd = ^Pz at both surfaces, and the induced charges 
cause a full depolarization field in the thin film, 5d ~ 
—4:TT(TindZ — —AttPzZ. On the other hand, if the ferro- 
electric thin film is placed between short-circuited perfect 
electrodes as depicted in Fig. [^b), the induced charges 
are fully canceled by free charges Cfree arising at both sur- 
faces of the electrodes, 6d — — 47r((7ind+Cfroo)2^ = 0. This 
geometric circumstance can be simulated with the dou- 
bly periodic supercell as depicted in Fig. djc), because 



the two electrodes act as two electrostatic mirrors fac- 
ing each other, and the mirrors make oppositely charged 
infinite mirror images beyond the electrodes. 

We can also introduce dead layers of thickness d be- 
tween the ferroelectric thin film and electrodes by con- 
straining the local soft-mode amplitudes to vanish {u = 
0) in these layers, as illustrated in Fig. E^d). With 
the dead layers, the infinite mirror images beyond the 
electrodes become j^ more sparse than images of the 
without-dead-layer configuration. Consequently, the free 
charges arising at the electrode surfaces decrease to 
Cfree = ^jTj<^ind, where I is the ferroelectric film thick- 
ness. This simulates short-circuited imperfect electrodes 
resulting in a depolarization field of 



Sd 



-An- 



l + d 



-,P.z . 



(26) 



We can also use a doubly periodic supercell with dead 
layers for this case. Physically, the depolarization field 
of Eg. (f26| can arise either from the presence of a 



dead layer in the ferroelectric near the interface, or 
from imperfect screening at the metal electrode, or 
both. We can define an cfFcctivc screening length for 
each of these effects, and we interpret the "dead-layer 
thickness" d of our model as corresponding to the sum 
of these two physical screening lengths. The screening 
length associated with the electrode interface appears 
in Eg. (f6) of Ref. 1\ and Eg. (1) of Ref. [S] an d is 
dis cussed for the SrRuOs/BaTiOs interface in Refs. 17|, 



18|, and 19|. Therefore, while the model does not 
explicitly incorporate information about the interface 
screening, this information is effectively included in the 
definition of the total screening length d in our model. 
Thus, for example, simulations at constant d for various 
film thicknesses can give the thickness dependence of the 
properties of capacitors with a certain interface structure. 
In the present MD simulations, the local soft-mode am- 
plitude vectors u in dead layers are fixed to zero by the 
infinitely large mass. This infinitely-large-mass trick is 
congenial to the Nose-Poincare thermostat for maintain- 
ing the Nose-Poincare Hamiltonian at zero. Moreover, 
this treatment also has another advantage in that the 
short-range interactions between the surfaces of ferro- 
electric thin film and the electrodes are automatically 
truncated. 

The depolarization field £d increases the total energy 
of the ferroelectric thin film hy —P ■ £d = iTtj^P^. To 
avoid forming a depolarization field in ferroelectric thin 
films, it is known that the films often develop striped do- 
main structures i^i2ii^t?3 ^^j^g introduction of the striped 

domain structure can eliminate some part of the energy 
increase An-j-^P^, because Pz becomes zero on average. 
However, the striped domain structure involves an energy 
cost in the short-range interaction V^^°'^*' {{u}) , because 
it has domain boundaries between which u has opposite 
direction zLz. The shorter the wavelength A of the striped 
domain structure, the weaker the depolarization field, 
but the higher the short-range interaction energy. The 
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FIG. 2: Schematic illustrations of ferroelectric thin films of 
thickness I unit cells (here 1 = 2). (a) Isolated thin film in 
vacuum, (b) Thin film sandwiched between short-circuited 
perfect electrodes. Doubly periodic boundary conditions for 
simulations of films sandwiched between perfect and imper- 
fect short-circuited electrodes are depicted in (c) and (d), re- 
spectively. Horizontal thick lines marked with "E" represent 
the electrostatic mirrors used to model electrodes. They are 
a distance d/2 away from the ferroelectric film surface (d = 
in (c), d = 1 in (d)). Each thin arrow represents a local dipole 
within a unit cell (a^ = 3.94 A^) of the BaTiOa crystal. Thick 
dashed lines enclose the periodic cell used for simulations. 



ground state of a ferroelectric thin film will be decided 
by a competition between the long-range dipole-dipolc in- 
teractions which favor a short-period domain structure, 
and domain-wall energy that arises from the short-range 
interactions and favors a uniformly polarized structure 
or a longer-period striped structure. In some previous 
^Qj.]jg^24j25|26 -j-j^g imperfect screening was mimicked with 
a parameter. On the other hand, our method with doubly 
periodic boundary condition does not reguire any param- 
eters, because the effect of imperfectness of electrodes is 
automatically and implicitly included in the long-range 
dipole-dipole interaction y^P'({it}). 



III. RESULTS AND DISCUSSION 

A. Bulk BaTiOs 

We first check the reliability of our MD program 
by comparing results of our simulations for bulk 
BaTiOa with earlier work based on the same effective 
Hamiltonianii^iii We used a system size of L^xLyXEz = 
16 X 16 X 16 and small temperature steps in heating- 
up (-1-5 K/step) and cooling-down (—5 K/step) simu- 
lations, with initial configuration generated randomly: 
(ua) = 0.07 A and (u^) _ (u^f = (0.02 A)^. We have 
also checked that there was no dependence of results of 
these simulations on initial configurations. The tempera- 
ture dependence of the homogeneous strain components 
(see Fig. ^, which are the secondary order parameters 
of ferroelectric phase transitions, exhibits the correct se- 
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FIG. 3: Average homogeneous strains e^x, Gyy, and e^z as 
a function of temperature in heating-up (+5 K/simulation, 
solid lines) and cooling-down (—5 K/simulation, dashed lines) 
simulations for a 16 x 16 x 16 supercell. Strains are measured 
relative to the LDA minimum-energy cubic structure with 
lattice constant 3.948 A. 



quence of phase transitions in BaTiOs known experimen- 
tally. Even under the negative pressure p = —5.0 GPa, 
the paraelectric to ferroelectric transition temperature 
Tc is underestimated at around 320 K in comparison with 
the experimental value of Tc — 408 K. Our estimates of 
Tc's agree fairly well with the ones reported in Ref. [ll| . 
The relatively weak first-order nature of the cubic-to- 
tetragonal phase transition in comparison with the first- 
order tetragonal-to-orthorhombic and orthorhombic-to- 
rhombohedral phase transitions is evident in the width 
of the temperature intervals of hysteresis (see Fig. [3]). We 
note that the ability to simulate time-dependent phenom- 
ena is one of advantages of MD simulations compared to 
Monte-Carlo simulations. 



B. BaTiOa ferroelectric thin- film capacitors 

We now simulate and analyze the behavior of epi- 
taxially grown films of BaTiOs on GdScOa sub- 
strates.^ In our simulations, we represent this with 1% 
in-plane biaxial compressive strain by maintaining the 
homogeneous strain r;i = 772 = —0.01 and rje = 0. In 
other words, we maintained the average lattice constants 



a and b at 0.99ao and angle 7 at 90°. We use super- 
cell sizes oi Lx X Ly X Lz = 32 x 32 x 2(Z -I- d) and 
40 X 40 X 2 (? + d) and simulate ferroelectric layers of thick- 
ness I sandwiched between two short-circuited electrodes 
with {d = 1) and without {d — 0) dead layers. This is 
accomplished through use of doubly periodic boundary 
conditions as explained earlier. 

Both the heating-up and cooling-down simulations are 



started with an initial configuration of {u^ 



0, 



(uz) = 0.07 A and {uD - (u„)2 = (0.02 A)^. In the 
cooling-down simulations, which start at a sufficiently 
high temperature, the initial configuration changes to an 



unpolarized one {{uz) = 0) during thermalization. We 
monitor the temperature dependence of (ua) and (u^) for 
thin films with thicknesses I = 15, 31, 127, and 255 with 
dead layers d = 1 and a thin film with thickness I — 32 
without dead layers {d — 0) (see Fig. |4] and animations 
in the EPAFS-^ ). The behavior of the film with no dead 
layer is the same in heating and cooling simulations. In 
contrast, for the films with a dead layer (d — 1), the tran- 
sition behavior exhibited by (uz) = is rather different 
in heating and cooling simulations, although the temper- 
ature dependence of (m^) is almost the same in both the 
kinds of simulations. In the heating-up simulations, the 
discontinuity in (uz) as a function of temperature marks a 
transition from a ferroelectric state with almost uniform 
out-of-plane polarization (Fig.[5lja)) to one with a striped 
domain structure (Fig. [SJb) and (c)). We find that this 
transition temperature, Ts{l, d = 1), exhibits a strong de- 
pendence on size I. We note that this transition is missing 
in the cooling-down simulations; just above Tg, striped 
domain structures appear and the stripes remain and be 
frozen at T < Tg. The temperature Tc{l, d = 1) at which 
•\/(uf)(T) exhibits a change in its slope marks another 
transition, namely from a striped domain phase to a para- 
electric phase. Tc{l,d = 1) depends relatively weakly 
on the film thickness. It should be mentioned that the 
results of heating-up simulations with a phase transition 
from the single domain state to the striped domain state 
at Tg below Tc agree with thcrmodynamical treatment 
of ferroelectric capacitors with dead layers by Chensky 
and Tarasenko. 



"^S" 



For films with d = 1, Tg is 150 and 210 K for / = 15 
and I = 31 respectively, which is lower than the bulk 
transition temperature (Tc ~ 320 K). However, for d = 1 
films with / = 127 and I — 255, Tg is enhanced to 520 and 
610 K respectively, well above the bulk Tc. In the infinite 
thickness limit (Z -^ 00), it appears that Ts{l,d = 1) 
tends to the Tc of thick films with no dead layer (d = 
0), since Tc is 650 K for / = 32 and d = 0. In the 
d — I cases with I = 127 and / = 255, the effect arising 
from the depolarization field weakens significantly, and 
the enhancement of Tg results from the in-plane biaxial 
compressive strain. In the d = case with I = 32, there is 
no depolarization field and enhancement of Tc by the in- 
plane biaxial compressive strain is effective even in very 
thin films. We note that y/{ul) and Vv"I) ^^^ distinct 
even at high temperatures (see Fig. 2^), indicating that 
the symmetry of the paraelectric phase is broken by the 
presence of the epitaxial constraint and the electrodes, 
as well as correlations between local dipoles and their 
images. 

For films with a dead layer {d = 1), the striped do- 
main structures appear in the cooling-down simulations 
at low temperatures for all values of thicknesses / ex- 
plored here (see Fig. [Hb) and (c) for the case oi I — 15 
with d — 1, and Fig. [6] for various I). As shown in Table [H 
the wavevector k of the striped domain, at which Uz{k) 
has the largest amplitude [^^(fc)!, exhibits an interesting 
dependence on thickness /. We have determined k for two 



TABLE I: Dependence of the wavevector k/2-K of the striped 
domain structure on thickness I in the thin-film BaTiOa ca- 
pacitor with a dead layer (d = 1) . 



/ 


32 X 32 X 2{l -f- d) 


40 X 40 X 2{l + d) 


7 


{ 4/32 3/32 } 


{ 5/40 5/40 } 


15 


{ 3/32 3/32 } 


{ 4/40 3/40 } 


31 


{ 2/32 2/32 } 


{ 3/40 2/40 } 


63 


{ 2/32 1/32 } 


{ 2/40 2/40 } 


127 


{ 1/32 1/32 } 


{ 1/40 1/40 } 


255 


{ 1/32 0/32 } 


{ 1/40 0/40 } 



supercell sizes, 32 x 32 x 2{l -f d) and 40 x 40 x 2{l + d), to 
identify supercell-size effects. It can be seen that, except 
for the data for I = 255, k tends to be along the in-plane 
{110} direction, consistent with earlier reports. ^'''■^^ The 
simulated striped domain structure for I — 255, which 
is parallel to the {100} direction, is likely to be an ar- 
tifact of the finite supercell-size: Lx x Ly = 32 x 32 or 
even 40 x 40 are too small to allow for the formation of a 
sufficiently thick {110} striped domain. The wavelength 
A = 27r/|fc| of dominant periodicity of the domain pattern 
is shown as a function of thickness / in Fig. [71 where it is 
evident that the thinner films have smaller A to avoid the 
stronger depolarization field, Eq. ((26)) . The fitting shown 
in Fig. [7] suggests a square-root dependence^Si on I (the 
result for I — 255 is not included in the fit). Extensive 
simulations at larger length scales would probably be re- 
quired to clarify further the dependence of the domain 
period of these striped structures on film thickness and 
dead-layer thickness. 

The stark difference in the behavior of (uz) in heating- 
up and cooling-down simulations hints that the (almost) 
uniformly polarized state and the (uz) = striped do- 
main states are frozen and thermal hopping between 
them may be almost impossible at low temperatures. To 
understand why both uniformly polarized and striped 
domain states are stable and thermal hopping between 
them are difficult, we investigated the effective potential- 
energy surfaces for striped domain structures of various 
stripe wavevectors k and various I for thin-film ferro- 
electric capacitors with and without the dead layer (see 
Fig. [5Ja)-(e)). Omitting surface relaxations in this anal- 
ysis may be reasonable because the surface relaxations 
are confined to the surface region of the ferroelectric thin 
films as shown in Fig. \Slc). It can be seen that thin- 
ner ferroelectric film have a shorter stripe wavelength 
A = 27r/|fc| in their ground states. As the thickness / is 
increased to I w 127, the ground state changes from the 
striped domain structure to the out-of-plane uniformly- 
polarized ferroelectric structure (fc = (000)). However, 
on the time scale of our simulations (w 1 ns), even at 
I « 255 there is no hopping from the striped domain 
metastable state to the uniformly polarized ground state 
(Fig. IH^d)). It can also be seen in Fig. [5Ka)-(e) that the 
magnitude of u^ which gives the minimum-energy ground 
state becomes larger, and the minimum energy gets 



deeper, as I increases, in good correspondence with the 
thickness dependence of Tq- The trend of k with I also 
shows good agreement with the simulated values shown 
in Table |T1 The simulated stability of the out-of-plane 
uniformly-polarized states against the energetically lower 
striped-domain states in thinner (/ < 127) films at low 
temperature seems to give support to the recent idea of 
elastic stabilization of a homogeneously polarized state 
in strained ultrathin films. •^" As shown in Fig. [IKa), the 
polarization switching in the epitaxially constrained film 
may be suppressed by the presence of a potential barrier 
that prevents hopping between the uniformly-polarized 
and striped-domain states. For I < 127 with d = 1, it 
is expected that a uniformly polarized film would evolve 
into a striped domain state, or vice versa, over a suffi- 
ciently long time at T < Tg. However, the time scale 
of the evolution might be very much longer than the 
present simulation time scale ('--^l ns). It might also be 
expected that, in the cooling-down simulations of films 
with I < 127 and d — I, the uniformly polarized state 
is obtained at T < T^. Instead, however, we find that 
stripes appear. A close inspection of the simulations 
shows that the stripes form slightly above Tg, initially 
in a somewhat disordered fashion, presumably because 
such a structure provides a good compromise between 
energetic and entropic considerations. The stripes then 
get frozen into place, and become better ordered, as the 
temperature is reduced below T < Tg. Conversely, in 
the case of d = (Fig. Hfe), depolarization field E^ — 0), 
the striped domain stricture does not appear during the 
heating-up and cooling-down simulations. This may be 
because, when £d = 0, there is no reason or chance to 
form a striped domain structure even just above Tq- At 
Tc, direct phase transition form paraelectric phase to 
uniformly polarized ferroelectric phase occurs. Then, be- 
low Tc, the system tends to be in its ground state, the 
uniformly polarized ferroelectric structure. 



C. Hysteresis loops 

A measurement of polarization typically involves use 
of a triangle-wave electric field for recording the ferro- 
electric hysteresis loops (inset of Fig. [T0| . The hystere- 
sis loops and coercive fields £c depend on the ampli- 
tude £q and frequency / of the applied fields. We simu- 
late hysteresis here using triangle- wave with steps (width 
A^ "-steps and height A£) as sketched schematically in 
Fig. [ini Thus, the frequency of the applied field in our 
simulations is / = A£/4 Aingtops^o- We used supercell 
sizes oi LxX LyX Lz ~ 16 x 16 x 2(Z + d) in simulations of 
hysteresis loops for ferroelectric thin-film capacitors with 
1% in-plane biaxial compressive strain and without con- 
straints of strain (namely, the "free" film) (see Fig. [TT|) . 
The temperature is maintained at 100 K through the 
simulations. For both the epitaxially constrained and 



"free" films, our simulations confirm that the imperfect 
screening of the electrodes decreases the coercive field as 
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the film thicknesses decreases, as described phenomeno- 
logically in Ref . @ . There is a large (order-of-magnitude) 
difference in the coercive field £c between the epitaxially 
constrained film and the "free" film. This may be because 
the compressive strain arising fi'om epitaxial constraints 
prevents the polarization switching, while the inclusion of 
inhomogeneous strain (i.e., acoustic displacements) eases 
the switching, as depicted in Fig. [HI The potential bar- 
riers themselves are lower in the "free" films than in the 
epitaxially constrained films (see Fig. [5]). We note that 
hysteresis loops for "free" film capacitors with I — 63 and 
I — 127 are very similar to the experimentally observed 
hysteresis loops of a ferroelectric capacitor with damaged 
electrodes that have "steps" and "plateaus" during po- 
larization switchings.^ This is because, in the "free" film 
capacitors with imperfect electrodes {d = 1), the con- 
figuration with out-of-plane polarization is no longer the 
ground state. In fact, the ground state has a nonzero 
in-plane polarization. Thus, the dipoles Z*u{R) have 
large in-plane components Z*Ux{R) and Z*Uy{R) in the 
hysteresis- loop simulations (and experiments), as evident 
in the snapshot shown in Fig. 1121 

Unfortunately, attempts to fit our results to the 
usually-assumed Kay-Dunn scaling of the coercive 
field £c with film thickness I of thicker films'*^ were 
unsuccessful, as were attempts to emulate the relatively 
weak dependence of £c on I for epitaxially grown 
high-quality ultrathin filmsj '^^i'^^ ''^ The experimentally 
observed values of coercive fields £c for ultrathin 
BaTiOa capacitors range from 200 to 500 kV/cm,^^'^^-^-^ 
while simulations of epitaxially constrained films 
largely overestimate £c, and those of "free" films 
slightly underestimate £c- This may be because the 
switching in real thin-film capacitors is a large-scale 
(> 100 nm) phenomenon involving defect-mediated 
nuclcation mainly at ferroelectrics-electrodes 
interfaces , ^^i'^^'^^ '^''-'^^-'^^ as well as the possibihty 
that the strain conditions may be intermediate between 
the cases of epitaxially constrained and "free" films. 
Such intermediate strain conditions may be achieved and 
will be simulated with MD in the future by introducing 
a mechanical boundary condition such as presented 
in Ref. 30|. In contrast to our case of ultrathin 
BaTiOa capacitors, it is well known that for ultrathin 
PbZra;Ti_.r03 (PZT) capacitors the coercive fields £c 
increase with decreasing film thickness I, and there is an 



argument whether this strong increase of £c is coming 
from compressive substrate-induced lattice strain"''* or 
not.**^ Constructing a first-principles Hamiltonian for 
PZT and simulations with this MD method will help us 
to understand this difference between BaTiOs and PZT. 



IV. SUMMARY 

We have developed a robust and highly efficient 
molecular-dynamics scheme, based on a first-principles 
effective Hamiltonian formulation, for simulating the 
behavior of the polarization in perovskite-type ferro- 
electrics. We have applied this approach to study 
BaTiOa ferroelectric thin-film capacitors, with special at- 
tention to the dependence on film thickness and choice 
of electric boundary conditions. We find that striped 
domain structures tend to form on cooling-down simula- 
tions when a ferroelectric dead layer is present near the 
electrodes, and we study the dependence of the domain 
period on the conditions of formation. We also study 
the hysteresis loops for capacitor structures, both with 
and without such dead layers, and we find dramatic dif- 
ferences in the hysteretic behavior for the cases of elas- 
tically constrained or "free" films. Our MD simulator 
f eram will be a powerful tool for further investigations 
of the physical properties of ferroelectric nanostructures 
that are relevant for a variety of potential device appli- 
cations. 
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FIG. 4: (mz) of heating-up (solid lines) and cooling-down 
(dashed lines) molecular-dynamics simulations of BaTiOa 
thin-film capacitors with short-circuited electrodes under 1% 
in-plane biaxial compressive strain for (a) thickness I = 
15 layer with dead layer d = \, (b) I = 31 with d — 1, (c) 
I = 127 with d = 1, (d) Z = 255 with d = 1, and (e) I = 32 
without dead layer (d = 0). Vv"!) are also plotted in (a)- 
(e). In (c)-(e), heating-up \/{ul) and cooling-down \/ {u1) 
\u%) is plotted only in (e), because 
and \/{Uy) are essentially identical 
in cases (a)-(e) (for both heating-up and cooling-down). Su- 
percells are of size 32 x 32 x 2(/ -f d). Animations of these 
cooling-down and heating-up simulations are also available in 
the EPAPS.-^ 
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FIG. 5: Snapshots at T = 100 K in heating-up ((a)) and 
cooling-down ((b) and (c)) simulations of ferroelectric thin- 
film capacitors of Z = 15 with d — 1. (a) and (b) are horizontal 
slices, (c) is a vertical cross section. Points of snapshots are 
indicated with "X" marks in Fig.|4ja). In horizontal slices, the 
-|-z-polarized and — 2-polarized sites are denoted by D and ■, 
respectively. In vertical cross sections, the dipole moments of 
each site are projected onto the a:«-plane and indicated with 
arrows. Layers which do not have arrows are dead layers. 
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FIG. 6: Horizontal slices of snapshots at 100 K in cooling- 
down simulations of ferroelectric thin-film capacitors with sin- 
gle dead layer (d = 1) of various thickness 1 = 7, 15, 31, 63, 
127, and 255. The -|-2-polarized and — z-polarized sites are 
denoted by D and ■, respectively. 
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FIG. 7: Calculated thickness I dependence of wavelength A 
of striped domain structures in thin film BaTiOs capacitors 
with a dead layer (d = 1). + marks are from 32 x 32 x 2{l + d) 
supercell calculations and x are those of 40x40 x2(Z+(i). Data 
oil <— 127 are fitted with A — cyl (dotted line). I = 255 data 
are omitted, because of their large supercell-size dependence. 
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FIG. 8: Effective potential surfaces of BaTiOs thin-filin ca- 
pacitors with short-circuited electrodes: (a)-(e), under 1% 
in-plane biaxial compressive strain arising from epitaxial con- 
straints; (f)-(j), without epitaxial constraints (i.e., for "free" 
films). The thicknesses of ferroelectric films and dead layers 
are indicated in each panel with I and d respectively. Total en- 
ergies as functions of Uz are compared among striped domain 
structures with wavevectors k parallel to (110). k = (000) 
corresponds to the uniformly polarized structure. The zero 
of the energy scale is placed at the total energy of the non- 
polarized Uz — structure. A negative pressure p = —5 GPa 
is applied to correct the underestimation in Tc ■ 
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(a) epitaxially constrained film (difficult to switcli) 




(b) "free" film (easy to switch) 



FIG. 9: Schematic comparison between the epitaxially con- 
strained film and the "free" film. In the epitaxially con- 
strained film, switching may have to climb over a potential 
barrier, but, in the "free" film, dipoles can be easily rotated 
and switching can go around a valley of the potential. 



E 




E 


Eo- 


1_ 


Afx«steps 




-''"^~ .f 


\,-'' \,' ' 






^ 


experiments (triangle wave) 
, f 






L 


"l — 1 simulations 


-Eo 






1 — 1 1 — 1 











FIG. 10: Schematic illustrations of triangle-wave electric field 
used to measure ferroelectric hysteresis loops experimentally 
(inset) and in the present simulations. 





60 
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(b) "free" film 
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FIG. 11: Calculated hysteresis loops for capacitors with (a) 
epitaxially constrained films, and (b) "free" films of vari- 
ous thickness I and with dead layer d. Supercell sizes were 
16x 16x 2(Z-|-d). Ai = 2 fs and nstcps = 50,000. For epitaxially 
constrained films. So — 4,000 kV/cm and A£ = 100 kV/cm 
are employed. For "free" films, £o = 400 kV/cm and 
A£ = 10 kV/cm are employed. 
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FIG. 12: Vertical cross section of a simulated ferroelectric 
"free" film capacitor with a single dead layer; 16 x 16 x (? = 
63, d= V). The snapshot was taken at the point marked "x" 
in Fig. Illf bl . The projection of the dipole moments onto the 
a;2-plane are indicated with arrows. 



